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Abstract. A new method for the formal solution of the 2D radiative transfer equation in axial symmetry in the presence of 
arbitrary velocity fields is presented. The combination of long and short characteristics methods is used to solve the radiative 
transfer equation. We include the velocity field in detail using the Local Lorentz Transformation. This allows us to obtain a 
significantly better description of the photospheric region, where the gradient of the global velocity is too small for the Sobolev 
approximation to be valid. Sample test calculations for the case of a stellar wind and a rotating atmosphere are presented. 
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1. Introduction 

Radiation is the major source of information about stars. A 
number of stellar properties can be obtained by comparing ob- 
served spectra with the synthetic ones calculated from model 
stellar atmospheres. A key problem in computing synthetic 
spectra is solving the radiative transfer equation in stellar at- 
mospheres, since the emergent radiation is formed in these re- 
gions. There exists a large number of methods for solving the 
transfer equation in a static one-dimensional case, which are, 
unfortunately, inappropriate for stars with stellar winds, rapidly 
rotating stars, accretion discs, and also for nebulae. Our main 
interest is to study various aspects of the stellar wind in hot 
stars, therefore we want to develop a method for solving the ra- 
diative transfer equation that is well suited for this case, where 
the symmetries enabling a one-dimensional approximation are 
broken. 

There are not many methods available to solve static multi- 
dimensional radiative transfer problems. For optically thin re- 
gions, the Monte Carlo method (Boisse 1990 1 may be used. 
On the other hand, in the optically thick regions it is possi- 
ble to solve the transfer equation using a method that employs 
the diffusion approximation (Kneer & Heasley[l979|. Neither 
of these methods are suitable for stellar atmospheres, where 
the transfer problem must be solved from optically thick re- 
gions to regions where the optical thickness is very small. The 
classical way of solving the radiative transfer equation in more 
dimensions is the long characteristics method (Cannon 1970 1. 
This method fully describes the radiation field, but the com- 
puter time necessary to obtain a solution may become long. For 



this reason, Kunasz & Auer ( 1988 1 developed the short charac- 
teristics method, which is the best currently available multi- 
dimensional method. There exist several applications of short 
characteristics methods in a Cartesian grid (Fabiani Bendicho 
2003 and references therein) and 2D axially symmetric geom- 
etry (Georgiev et. al l2003l Dullemond & Turolla d2000l de- 
veloped a "rotating plane" method for axially symmetric prob- 
lems based on the short characteristics method, too. An effi- 
cient approach is to use adaptive grids following Folini et al. 
( 2003 ) and Steinacker et al. J2003i . A review by Auer (|2003b I 
provides an excellent insight into the grounds for astrophysical 
multidimensional radiative transfer. 

Another possibility is to apply the finite element method, 
which has been recently used for multidimensional radiative 
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transfer by, e.g., Richling et al. (2001 1. However, the finite ele- 
ment method is not used very often for radiative transfer in stel- 
lar atmosphere studies due to its convergence problems caused 
by the about ten orders of magnitude changes of physical quan- 
tities (e.g. opacity, density) in stellar atmospheres. Dykema et 
al. 1^996) tried to overcome the instability and convergence 
problems using a modification of the finite element method, 
namely, the discontinuous finite element method. The discon- 
tinuous finite element method was also used in our previous 
paper (Korcakova & Kubat 2003 hereafter Paper I) for one- 
dimensional radiative transfer. 

If the velocity field is present, the situation is more com- 
plicated. The changes in opacity and emissivity along a ray in 
spectral lines can be large due to the Doppler shift, and we have 
to take this into consideration. In the continuum, we can simply 
use the static equation, since the Doppler shift has only a small 
effect on the opacity and emissivity coefficients. For simplic- 
ity the methods assuming monotonic velocity field have been 
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used very often. The multidimensional problem with velocity 
fields can be solved in the observer frame (the frame connected 
with the center of a star) or in the comoving frame (the frame 
connected with the outflowing material). Radiative transfer in 
more dimensions has usually been solved in the observer frame 
(e.g. Carlsson & Stein 2003). One of the firsts works to do this 
is Mihalas, Auer & Mihalas dl978> . who were able to solve 
a periodic velocity field in two-dimensional planar geometry. 
Recently, van Noort et al. (|2002 2003 1 solved this problem in 
two dimensions. Their code is able to cope with both spherical 
and cylindrical geometry for Cartesian coordinate systems. The 
observer frame is appropriate only for small velocity gradients. 
The velocity difference between neighbouring depth points has 
to be smaller than several times the thermal velocity. If the ve- 
locity gradient is large, the Sobolev approximation (Sobolev 
I1947> is often used. This approximation was used for a 3D case 
by Folini et al. J2003> . Both cases of small and large velocity 
gradients can be solved in the comoving frame. For a 3D ac- 
cretion disc, the radiative transfer in the comoving frame was 
solved by Papkalla JT9951 . 

However, there exists another effect that is usually ne- 
glected or taken into account in an approximate manner, 
namely the rotation of stars. Stellar rotation is usually taken 
into account as a convolution of the line profile, obtained from 
the plane-parallel static solution, and of the rotation profile (cf. 
Gray 119761 . In this technique one must make use of analyti- 
cal expressions for limb darkening, which do not give a correct 
description of the angular dependence of the specific intensity. 
It also does not take into account the fact that limb darken- 
ing is strongly frequency dependent across the line profile (cf. 
Hadrava & Kubat 2003 ). In order to describe the radiation field 
and its influence on the stellar atmosphere correctly, it is also 
necessary to take into account the Doppler shift of lines in a ro- 
tating atmosphere for the transfer of radiation along nonradial 
rays. 

In this paper we present a new method to formally solve 
the radiative transfer equation in axial symmetry, which does 
not employ the Sobolev approximation. Our method is based 
on the Local Lorentz Transformation method (LLT) described 
in Paper I for the one-dimensional case. For studies of stellar 
wind, accretion discs or stellar rotation we need a method that 
uses a more general geometry than plane-parallel or spherical. 
However, it is not necessary to treat the whole three dimen- 
sional space in detail. We employ axial symmetry. The results 
will be more accurate than in a plane-parallel or spherical ge- 
ometry and the calculation will be faster than for the full 3D 
problem. 

In the first part of this paper we present our axially symmet- 
rical code. The geometry of the model is described in detail, 
followed by some tests of our code. We compare these results 
with a spherically symmetric model atmosphere from Kubat 
(2003 1. We also present some test calculations with a velocity 
field. Our line profiles are compared with convolved profiles 
for the case of stellar rotation. Some results for the stellar wind 
are also shown. Then we test the dependence of the computing 
time on the number of grid points. In the last part, we comment 
on the possibilities of the application for our method. 



2. Method 

We solve the radiative transfer equation for the specific inten- 
sity of radiation / 

n ■ VI (r, n, v) = ~x (r, v) [I (r, n,v) - S (r,v)] , ( 1 ) 

where x is the opacity, S is the source function, n is the direc- 
tion of radiation, r is the radius vector, and v denotes the fre- 
quency. We assume axial symmetry and we allow for a nonzero 
velocity field. 

The basic idea of this method is to solve the radiative trans- 
fer problem not in the whole star, but in separated planes inter- 
secting the star. 

2.1. The spatial grid 

Let us consider the spherical coordinate system (r, 8, <p). We 
choose as the axis of symmetry 8 = (see Fig.[0. First, we 
introduce the discretization of the radial distance r and angle 8. 
Due to the axial symmetry, physical quantities do not depend 
on the angle <p. The grid is chosen to give the best description 
of the system, depending on the studied object. As an exam- 
ple, for the study of stellar winds, the grid of angles 8 can be 
equidistant. On the other hand, for accretion discs it must be 
finer near the equatorial plane, where the disc resides. In radial 
distance r we choose the grid (very similarly to a ID problem) 
to be equidistant in the logarithm of the radial optical depth. 
A suitable choice is about 5 points per decade. We assume the 
opacity and emissivity of the stellar material to be known at the 
grid points. 

To reduce the 3D problem to a 2D one, we do not solve 
the radiative transfer equation directly in the 3D grid, but in a 
set of "longitudinal planes" intersecting the star parallel to the 
plane <p — (we "slice the star" - see Fig.Q. In thin stellar at- 
mospheres it is favourable to choose a finer grid of longitudinal 
planes closer to stellar limb, since the limb darkening is better 
described by such a choice. 

The primary grid described above helps us only to inter- 
polate the quantities to the selected planes and finally to cal- 
culate the emergent flux. To solve the transfer equation we 
need a slightly modified grid. For each longitudinal plane we 
choose the polar coordinate system and define a grid of con- 
centric grid circles and radial grid lines (see Fig.^. The grid 
of concentric circles corresponds to the original 3D grid in the 
planes, which intersect the opaque stellar core where we do 
not solve the transfer equation and where we apply the lower 
boundary condition following from the diffusion approxima- 
tion. However, choosing an appropriate grid for planes which 
do not intersect the opaque stellar core may become difficult, 
since the grid chosen must correspond to the geometry of the 
problem as well as resolving the velocity field. This is shown 
in Fig.|2l We must provide an appropriate grid for both geomet- 
rically thin and extended atmospheres. 

For geometrically thin atmospheres the grid of radial dis- 
tances in the longitudinal plane is chosen using intersection 
points with grid circles in the primary 3D grid (left panel in 
Fig-GJ- F°i" an extended atmosphere it is better to define grid 
points where the plane intersects the radial grid lines (right 
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Fig. 1. Coordinate system and the set of longitudinal planes. 
These planes are parallel to the plane = 0. Note that some 
of them intersect the opaque stellar core, where the diffusion 
approximation is valid, while the rest does not. 

panel in Fig.13. The reason why we need an appropriate grid 
in the planes that do not intersect the inner boundary is the ne- 
cessity to ensure sufficient spatial resolution for the radiative 
transfer equation solution. This may be difficult in the central 
region of the plane. Moreover, similar problem involves by the 
velocity field. The velocity gradient in the planes far from the 
center can be large, so it is necessary to have a sufficient num- 
ber of grid points. 




Fig. 2. Choosing the grid in longitudinal planes that do not 
intersect the opaque stellar core. For geometrically thin (left 
panel) atmospheres, the grid points are defined using intersec- 
tions of a longitudinal plane with concentric circles. For ex- 
tended atmospheres (right panel), intersections of the longitu- 
dinal plane with radial lines are used for the definition of grid 
points. 



The opacity, emissivity and the source function are inter- 
polated to the new coordinate system. Linear interpolation is 
used here, since the grid is not orthogonal. The radiative trans- 
fer equation is solved for each longitudinal plane independently 
(Section l2~3l . The whole radiative field is obtained by rotating 
the separate planes around the axis of symmetry (Section lZ4l . 

2.2. The frequency grid 

The frequency grid is to be chosen to enable the most efficient 
description of the radiation field. We set the frequency interval 
using both the largest line width and the total Doppler shift 
caused by the global motion. For example, for Doppler profile 
this interval is 



v € <v (1 - (Voo + v ro , )/c) - 2Av D , 

V (1 + (Voo + v rof0 )/c) + 2Av D >, (2) 



where Voo is the terminal velocity, v rolo the rotation ve- 
locity in photosphere. Avp denotes the broadest Doppler 
halfwidth, which corresponds to the highest temperature. For 
very extended atmospheres or planetary nebulae (optically thin 
medium) we must take into account that the regions with — v M 
and v M "see" each other. In that case, we must multiply Voo 
in Eq. (|2ji by two. We take the frequency step to be equidis- 
tant for the case of the Doppler line profile. This step must be 
small enough to describe the line profile variations in the region 
with low temperature. The step is determined in the following 
way. We first typically set 5 points per line equidistantly at the 
depth point with the narrowest line profile (which is usually the 
depth point with the lowest temperature). Then we cover the 
whole frequency interval (|2j with frequency points using this 
frequency step. This frequency grid is then used in other depth 
points where the line profile is broader. For a Voigt profile we 
need to extend the frequency interval (|2ji, but the frequency step 
in the extended part of the interval may be larger. 

2.3. Solution in longitudinal planes 

In the given longitudinal plane we introduce the polar coordi- 
nate system using grid circles and radial grid lines defined in 
the preceding subsection ( 12. U . The values of temperature, den- 
sity and velocity are known at the grid points. First, we solve 
the transfer equation in the plane starting at the outer bound- 
ary. Once the radiation field in the direction towards the cen- 
ter is known, we can continue to solve the transfer equation in 
the opposite direction, i.e. from the inner boundary towards the 
outer one. At the inner boundary the diffusion approximation 
may be taken for planes intersecting the opaque stellar center as 
a boundary condition. In other planes (that do not intersect the 
stellar center), we adopt as the lower boundary condition the in- 
tensity taken from the previously calculated solution from the 
outer boundary inwards. 
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2.3.1 . The solution from the outer boundary to the 
central regions 

We begin the calculations at the outer boundary (stellar sur- 
face), where the boundary condition (i.e. the incoming intensity 
/) is known. At each inner grid point we choose several rays per 
quadrant (see Fig. [3}- The rays start at the outer grid circle and 
end at a given grid point. Along these rays we solve the transfer 




circle 



Fig. 3. Diagram for the solution of the radiative transfer equa- 
tion in the longitudinal plane starting at the upper boundary to 
the central regions. The rays along which we perform the for- 
mal solution are plotted using the dashed line. 

equation. The angle distribution of these rays may be the same 
as often used in the case of the plane-parallel geometry, where 
the angle cosines fi = cos a (a is the angle between the ray and 
the radial direction, see Fig. [3ji are chosen to be the roots of 
Legendre polynomials in the interval (0, = 0.8872983346, 
0.5, 0.1 127016654) to ensure better numerical accuracy of the 
angle integration (Press et. al ll986l section 4.5). Three rays at 
a given point is usually sufficient, because the whole radiation 
field is obtained by summing the information from all longitu- 
dinal planes intersecting the given point (see Fig.0. However, 
in some situations it is better to use more (up to 9) rays per 
quadrant to overcome a numerical error introduced by the nec- 
essary interpolation of intensity described below. As one can 
see in Fig.|3] the rays start at the preceding grid circle. This 
means that the rays do not usually start at a grid point and that 
they may also intersect some radial grid lines. 

The solution diagram is illustrated in Fig.|4] We perform a 
linear interpolation of the source function and opacity to obtain 
their values at points A, B, and C and of the incoming intensity 
for the value at point A. The optical depth difference At is cal- 
culated along the ray between the individual intersection points 
(abscissas AB, BC, and CD) using the linear approximation, 

At ( ab) = XA + 2 XB ks {AB) , (3) 

and similarly for At(bq an d ^ t (cd>- Here;^ and^g are opac- 
ities at respective points and As(^b) is the geometrical distance 
between points A and B. 




Fig. 4. The scheme for the integration along the ray for a de- 
termination of the specific intensity at the point D. The spe- 
cific intensity at the grid points of outer grid circle is known. 
The value of the specific intensity at the point A is determined 
by linear interpolation from values at neighbouring grid points. 
The values of the source function and opacity in the points B 
and C are obtained by linear interpolation. 

We solve the equation of radiative transfer by parts between 
all intersection points along the ray. We assume a linear depen- 
dence of the source function on the optical depth between the 
intersection points. For the interval AB the solution is 

S(t)e [ - {AT ^-' )] dt, (4) 

o 

and similarly for intervals BC and CD. The final intensity 1(D) 
at the grid point D is thus determined by three successive ap- 
plications of the equation 0. In this manner we obtain the spe- 
cific intensity in the downward solution at every grid point. It 
may, of course, happen that the geometric distance along the 
ray in the given cell is very small (the ray near the point B in 
the Fig. 0I. In this case we do not solve the transfer problem in 
this cell and we simply set 1(g) = 1(A)- Doing this, we eliminate 
a numerical instability. Usually, it is enough if the condition 
| cos(7r/2 - 9 — a)\ > 10~ 7 (a has the same meaning as in Fig. [5} 
is fulfilled, even if it is much more accurate to estimate it using 
the optical depth difference. 

There are several reasons why we prefer only a linear de- 
pendence of the source function on other higher order approx- 
imations. First, it is numerically stable and second, it is much 
faster. The higher order optical depth interpolation may some- 
times lead to numerical errors by adding new extrema (cf. 
Auer 2003a I. They may become significant especially in mov- 
ing medium. There is only one reason why linear interpola- 
tion should not be used: because the diffusion approximation 
in deep optically thick layers would be inaccurate. Since we do 
not solve the transfer problem using the short characteristics 
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method, this inaccuracy is not large in our case. Our character- 
istics are longer and the transfer equation is solved in several 
steps, so the information from the farther regions is naturally 
included. Therefore we chose the safer linear interpolation. 

2.3.2. The solution from the central regions to the 
outer boundary 

The upward solution is very similar to the previous step. The 
procedure is depicted in Fig.|5] We lead the rays to the grid 
points using the same angles as in the previous case. The rays 




Fig. 5. Diagram for the solution of the radiative transfer equa- 
tion from the center of the star to the outer boundary. The rays 
along which we solve the formal solution are plotted using the 
dashed line. 

start either at the preceding grid circle (closer to the center) or 
at the same circle as the grid point. At the intersection points we 
interpolate the opacity and the source function. Between these 
points we calculate the optical depth using (|3} and solve the 
transfer equation using Q as in the previous case of downward 
integration. 

From Fig.^ one can see that there exist planes that do not 
intersect the inner boundary region. For these planes we adopt 
the intensity calculated from the previous step (solution from 
the upper boundary to the stellar center) as the lower bound- 
ary condition. However, the solution of the transfer equation in 
the central grid circle must be performed with care. The situa- 
tion is shown in the Fig. [6] First, we solve the radiative trans- 
fer equation from the intersection point A to the center of the 
abscissa AC (point B) and then to point C. The physical quan- 
tities at point B we obtain by linear interpolation from points 
A and C. The optical depth difference is calculated using 
and the specific intensity is determined using (0}. The value of 
the intensity at point C determines the lower boundary condi- 



tion for further solutions outwards. For planes, which intersect 



c 




Fig. 6. Diagram for a solution of the radiative transfer equation 
in the central region in the planes, which do not intersect the 
region of validity of the diffusion lower boundary condition. 

the region of validity of the diffusion approximation (the stel- 
lar core), the situation is easier. We simply use the appropriate 
lower boundary condition. 

2.4. The full radiation field 

To obtain the whole radiation field we take the advantage of 
the symmetry of the problem. We know the radiation field at 
grid points in all directions lying inside the longitudinal planes. 
The radiation field in the whole star is then obtained by rotat- 
ing the longitudinal planes around the axis of symmetry 6 = 
(see Fig.0. This gives a sufficient description of the specific 
intensity. The emergent radiation flux towards the observer is 
calculated by integrating specific emergent intensity over the 
stellar disc. 

2.5. Velocity field 

The radiative transfer equation in moving media can be solved 
either in the observer frame or in the comoving frame. Since 
the coefficients of emissivity and opacity are angle dependent 
in the observer frame, we solve this problem rather in the co- 
moving frame. This frame, which is coupled with the moving 
medium, is generally a non-inertial frame. The form of the ra- 
diative transfer equation in the comoving frame we can obtain 
either from the general relativity or from the special relativity 
by assuming a set of Local Lorentz Transformations. 

In the Paper I we introduced an LLT (Local Lorentz 
Transformation) method for a solution of the radiative trans- 
fer equation in one-dimensional moving media based on an ap- 
plication of Lorentz transformations on cell boundaries and a 
static radiative transfer equation solution between them. The 
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Fig. 7. Diagram describing the determination of the whole ra- 
diation field. The rotation axis for longitudinal planes is per- 
pendicular to the page. Longitudinal planes (projected as lines 
in this figure) are rotated along the axis to intersect at a point R. 
Due to the axial symmetry, the radiation field for all directions 
is obtained by a rotation of all longitudinal planes around the 
axis of symmetry. 

method was tested by a comparison with an exact solution of 
the radiative transfer equation using the discontinuous finite el- 
ement method (see Fig. 8 in Paper I), which solves the comov- 
ing transfer equation including the d/dv (Doppler) and d/dfi 
(aberration) terms. Here we apply the same idea for the case of 
a two-dimensional arbitrarily moving medium. 

In our method, we replace the solution of the radiative 
transfer equation in the comoving frame by a set of Local 
Lorentz Transformations and a solution of the transfer equation 
in corresponding local inertial frames. In these inertial frames 
we solve a static radiative transfer equation, since it is Lorentz 
invariant. To do this, we must consider the constant velocity 
of every cell and the change of the velocity we allow only 
at the cell boundary (see Fig. [8jl. Two other conditions must 
be fulfilled there. First, velocity with respect to the observer 
frame must be low to neglect aberration. This permits us to 
solve the radiative transfer along straight characteristic lines. 
This assumption is valid, if v/c <K 0.01 (Mihalas et al. 119761 
Hauschildt et al. 19951. Second, let us assume a steady-state 
fluid. Due to this assumption we can neglect the time delay be- 
tween different parts of the medium. The correctness of this 
approach has been discussed in more detail in Paper I, where 
this method was tested for a plane-parallel geometry. 

The basic solution scheme within cells, which was de- 
scribed in Section l2~3l is not affected by the presence of the 
velocity field. We project the velocity field to the longitudinal 
planes. At cell boundaries we perform the Lorentz transforma- 




Fig. 8. The scheme for the inclusion of the velocity field. 
Velocity within each cell is constant and a change of velocity 
is possible only at the cell boundaries. 

tion of frequency v (see also Eq. (24) in paper I) 

(Av is the velocity difference between cells). Since the velocity 
field is not relativistic, we can simply use the classical Doppler 
law, which speeds up the solution. The Lorentz transformation 
of intensity I' (v') = I(v) (v' 3 /v 3 ) (see Eq, (23) in Paper I) at 
the cell boundary has a negligible effect, and we do not take 
it into account, since the intensity is proportional to the third 
power of the frequency ratio. This is close to one in most stel- 
lar applications. Since the equation of the radiative transfer is 
Lorentz invariant, we solve the static equation of the radiative 
transfer @ within each cell. 

Since we solve the transfer equation in the comoving frame, 
we have to do one additional step. We know the source func- 
tion and opacity in the frequency grid points in the rest frame 
coupled with a given cell. The incoming intensity from a 
neighbouring cell is known at different frequencies, which are 
Doppler shifted according to the velocity difference between 
the cells. This means that we must interpolate the intensity to 
frequency points coupled to the rest frame of a cell. 

To ensure a correct treatment of line transfer we have to 
guarantee that the frequency shift at the cell boundaries due to 
Doppler effect is smaller than a quarter of a Doppler halfwidth 
(see Paper I). If it is not, we must make the grid finer. This 
condition must be fulfilled even if we calculate with the Voigt 
profile since the center of the line has a Doppler core. This is 
the only limitation of this method. If the velocity gradient is 
too high, we must refine the grid, which slows down the calcu- 
lation. 

3. Test calculations 

Tests of the method are based on a model atmosphere of a main 
sequence B star with effective temperature = 17 x 10 3 K, 
gravitation acceleration logg = 4.12, and radius 3.26/?©. We 
obtain the state parameters, electron density and temperature, 
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using the hydrostatic spherically symmetric model atmosphere 
code ATA (Kubat 2003 1. We consider no incoming radiation as 
the outer boundary condition and the diffusion approximation 
as a lower boundary condition. 

3. 1. Static case 

For a basic comparison we took a model atmosphere calcu- 
lated using the static computer code ATA. The radiative trans- 
fer in ATA is solved using the long characteristics method and 
Feautrier variables (see Kubat ll993l[T994i r2003 1. We compare 
the result from the ATA code with the flux obtained from our 
new code (Fig.|9jl- The line profile was chosen to be Doppler 
here and the input parameters in the axially symmetric code 
are independent of 9. This comparison for this simplest case 
proves the basic correctness of the new code. 




5e-05 1 1 1 

4.56e+14 4.565e+14 4.57e+14 

v 



Fig. 9. The comparison of Eddington flux in the Ha line from 
the static spherically symmetric code ATA (crosses) and the 
axially symmetric one (solid line). 



limb darkening 




Fig. 10. Theoretical limb darkening obtained using our code 
across a half of the stellar disk in the Ha line. R+ denotes the 
stellar radius. 







limb darkening in continuum 
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Fig. 11. Limb darkening in the continuum. Crosses denote data 
computed from our model. The fitting functions f(x) and g(x) 
which express the commonly used limb darkening laws are de- 
fined by equations and 0. 



Limb darkening is also often described by the law (Gray 1 976i . 



3.1 .1 . Limb darkening 

An important result arising from the solution of the transfer 
equation is the limb darkening law. In Fig.^|we plotted the 
limb darkening law across the Ha line profile calculated from 
our code. We choose as an x-axis the distance from the cen- 
ter of star in star radius units (x = 1 for r = where /?* 
is the stellar radius). At this scale the effect of limb darken- 
ing is more clearly visible. To show this effect in more detail, 
we extracted the dependence of the specific intensity on the 
distance from the center of the stellar disc for a continuum fre- 
quency (Fig.ll 1> and for the central frequency of the Ha line 
(Fig.ll2>. The obtained data for the continuum are compared 
with functions usually used to express the limb darkening law. 
The first function (f{x)) is adopted from Allen's Astrophysical 
Quantities dl963> . 



I(x) = l-a-b + a{\- x z ) 1/z + b(l - x l ) = f(x) 



(6) 



Using the least squares method we obtain the parameters a = 
0.55 ± 0.02 and b = -0.20 + 0.01 to fit our results to Eq. 



I(x) = (l-e) + e(l-x 2 y /2 = g(x) 



(7) 



Fitting this case we obtain the value of the parameter e = 
0.277 ± 0.008. 

Figure El snows the limb darkening in the center of the 
Ha line. As one can see, brightening instead of darkening is 
observed near the stellar limb, which corresponds to the effect 
of flash spectra in solar chromosphere. The same result was 



obtained also by Hadrava and Kubat ( 2003 ). 

3.2. Stellar wind 

The case of the moving atmosphere is first tested for a spheri- 
cally symmetric stellar wind, for which we adopt the beta law 
for the dependence of the wind velocity v on radius r (see, e.g., 
Lamers & Cassinelli 1999b 



v(r) = Vo<, < 1 - 



1-1^ 



(8) 
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limb darkening in line 



B 1 -005 



x[R ,] 

Fig. 12. Center-to-limb variations at the central frequency of 
the Ha line. Note the apparent weak limb brightening before 
final darkening at the stellar limb. 

We choose a velocity in the photosphere of Vr = 200km ■ s _1 
and a terminal velocity Voo = 2000 km ■ s _1 . 

To check the consistency of our new method we compare it 
with the plane-parallel method described in Paper I. The plane- 
parallel radiative transfer equation is solved using the discon- 
tinuous finite element method. The velocity field is included 
using the Local Lorentz Transformation similarly to this paper. 
In Paper I the Local Lorentz Transformation method was com- 
pared with the solution of the radiative transfer equation using 
the discontinuous finite element method, where the frequency 
term was consistently included into the solution as an indepen- 
dent variable. Since the velocity field is non-relativistic, calcu- 
lations with aberration give the same results as without it (in 
agreement with Mihalas et al. 119761 Hauschildt et al. I1995> . 
The result from plane-parallel geometry is plotted using the 



0.95 




4.565e+14 



1 D LLT 

1 D LLT, without rim intensity ■ 
4.57e+14 



4.575e+14 



Fig. 13. The comparison of the Ha line profile obtained from 
our axially symmetric model (solid line) and a plane-parallel 
model (denoted by ID LLT). The dotted line is obtained using 
the same model, but the radiation with fi < 0.26 is not taken 
into account to approximate the curvature of the atmospheric 
layers (see also Fig.ll4>. The beta velocity law parameter/? = 1 
for all profiles. 

dashed line in Fig.^] As we can see, the agreement with the 
new (axially symmetric - solid line) results is not good. The 
reason is in the geometry used. Rays with large angle (low /i) 
in a curved atmosphere (spherically symmetric or axially sym- 
metric) do not even touch regions that are optically thick in 



continuum, whereas for a plane parallel geometry all rays with 
ju > finally reach the continuum optically thick regions (see 
Fig. 1141 . As a consequence, there is more continuum radiation 
in a plane parallel geometry, which results in shallower line 
profiles. This has been well known for a long time (see, e.g., 
Kunasz et al. 1975 1. This is why we plot another line profile, 



spherical symmetry 



plane-parallel 



Fig. 14. Intensity obtained from plane-parallel and axially 
symmetric geometries differ especially for very sideways (al- 
most tangent) rays. 



which also is obtained using the plane-parallel geometry, but 
in which the intensity incoming from directions with // < 0.26 
is set to zero. We can see a very good agreement in the blue 
wing of the line profile, where we see central parts of the star 
and where the assumption of plane-parallel atmosphere is ac- 
ceptable. However, the sphericity effects, limb darkening and 
limb brightening influence the red wing of the line. A different 
profile from a axially symmetric code is in agreement with the 
fact the line profile obtained in spherical geometry is different 
from that obtained in a plane parallel geometry. The difference 
depends on the sphericity effects on the temperature structure 
(see, e.g., Kunasz et al. 1975, Gruschinske & Kudritzki 1979, 
Kubat 1995, 1999). 

In Fig.^]we compare the line profiles for three different 
values of the parameter (3 (0.5, 1 and 2). To show the possibil- 
ity of our code to solve a more general velocity field we plot 
in Fig.^]the line profile obtained using a decelerating veloc- 
ity field. We choose a linear dependence of the velocity on the 
radial distance in a logarithmic scale. The photospheric and ter- 



s~ l and 200 km 



s respec- 



minal velocities are 2000 km 
tively. 

The most important point is the velocity gradient in the 
region of line formation. For this reason, lines corresponding 
to different values of the /3 (in Fig. 1151 parameter have dif- 
ferent profiles, even if the values of the photospheric and in- 
finity velocity are the same. We do not see the P Cygni pro- 
file, since the input parameters are based on the hydrostatic 
model, which produces a geometrically thin atmosphere. To 
obtain a P Cygni profile, not only enough high velocity gra- 
dient, but a sufficiently extended line formation region are nec- 
essary. Consequently, only a blue shifted absorption profile ap- 
pears in our case. Note that similar line profile shapes were 
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rotation profile G(v) is equal to 




P = 0.5 
p = 1 
P = 2 
deceleration 



Fig. 15. The profile of the Ha line for the case of a stellar wind 
for three values of parameter f3 = 0.5, 1,2 (see Eq. (jSJl). For 
comparison, the line profile affected by decelerating velocity 
field is also plotted (fine dotted line). 



obtained by Noerdlinger & Rybicki ( 1974 1 for a plane-parallel 
expanding atmosphere. 

3.3. Stellar rotation 

To test the stellar rotation we assume a spherical star and we 
consider the power-law dependence of the rotation velocity on 
the radius, 



v(r) = v(/?«) 



R, 



(9) 



where R* is the stellar radius and v(/?,) is the rotation velocity 
in the photosphere. We choose the parameter j to be equal to 1, 
which corresponds to equatorial discs formed by a stellar wind 



(Kroll & Hanuschik 1997) using conservation of the angular 
momentum. Stars which rotate near their critical velocity are 
far from being spherically symmetric. For example, the ratio 
of the polar and equatorial radii of a Eridani was determined 
to be about one half (Domiciano de Souza 2003 1 or even 1/5 
(Jackson et al. 2004 1. In these stars the approximation of spher- 
ical symmetry fails and we must solve the transfer problem in 
a more general geometry. However, although gravity darkening 
is neglected in our test case, our code is able to easily handle 
this effect. 

We compare line profiles calculated using our model with 
a rotational velocity field and detailed radiative transfer with 
other possibilities of calculation of the rotationally broadened 
profile. 

The standard and most commonly used way is to solve a 
detailed static plane-parallel radiative transfer equation for a 
static atmosphere and to convolve the resulting profile with a 
rotation profile using the relation (see Grav ll976l . 



= H(v) * G(v) 



f 



H(v - Av) G(Av) dAv. 



(10) 



Here, F v is the flux for a given frequency, F c is the continuum 
flux, H(v) is the normalized flux from the nonrotating star. The 



2(1 - e) Vl-(^) 2 + \ne ^-{^f 
G(v) = - v w v w (11) 



Av„ 



K *(l-f) 



in the interval |Av| < Av max , and G(v) - elsewhere. This 
expression assumes the limb darkening law in the form and 
the parameter e is the same as in the Eq. @. 

A more exact and also computationally more expensive ap- 
proach also uses the emergent radiation from the static atmo- 
sphere, but takes into account the angle dependence of the spe- 
cific intensity. The resulting profile is then calculated by inte- 
grating the specific intensity across the stellar disc. This ap- 
proach will be called "integrated static profile" hereafter. In 
these two latter approaches the radius dependence of the ro- 
tational velocity cannot be taken into account and the pho- 
tosphere is tacitly assumed to rotate as a rigid body. 

The most exact solution is to calculate the emergent ra- 
diation using a full solution of the radiative transfer equa- 
tion in a moving atmosphere, as has been done using our 
method. Results for the rotation velocity described by Eq. 
with v(R*) = 0.2v c are shown in Fig.^] Here, v r = -\/GM7jR^ 
is the critical rotational velocity, which is equal to 540 km ■ s _1 
for our case. In this figure we plotted the profile of the Ha line 




static case 
integrated static profile 
with velocity field 
convolution e=0 
convolution e=1 



0.82 - 




Fig. 16. The comparison of the static Ha- line profile and ro- 
tation profiles obtained using different methods for a rotation 
velocity 0.2 of the critical rotation velocity. The line profile cal- 
culated exactly using our 2D radiative transfer code with a cor- 
responding velocity field (|9} is plotted using the full line. The 
line profile that takes into account the surface velocity field, 
but with the local line profiles those calculated from the solu- 
tion of the static atmosphere, is plotted using the dashed line 
(integrated static profile in the figure). Line profiles calculated 
using the convolution of a static profile and the rotation pro- 
file with e = and e = 1 are plotted using the dotted line and 
dashed-dotted line, respectively. In addition, the profile calcu- 
lated for the static case (no rotation) is plotted using the thick 
dotted line. 

in the static case, the profile obtained from convolving the static 
profile with the rotation profile JTTt using values e = and 1, 
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the one from the integration of the static profile over the stellar 
disc, and the "true" profile calculated with the full influence of 
the rotation velocity field. The ratio of the equivalent width in 
the static case and the flux obtained from integrating of the line 
profile over the disc is of the order of 10~ 5 , which represents 
good numerical accuracy. 

Both line profiles obtained using convolution dlOi are 
deeper than those calculated from our code. The difference be- 
tween line depth obtained by convolution and line depth ob- 
tained from our code (7 oul - - 1)1 {I - I om ) is about 5% for pa- 
rameter e — and almost 20% for e = 1. The error of five 
percent is not too large, since sometimes the observed spec- 
tra have a lower accuracy. However, for high S/N spectra one 
may have an accuracy better than 1 %, which makes the error 
of 5% significant. There is no doubt about the detectability of 
a 20% difference. The commonly used value of the parameter, 
e ~ 0.6, yields remarkable differences as well. This difference 
is a consequence of the dependence of limb darkening on fre- 
quency across the line profile. As we can see from Fig.^ll the 
continuum intensity decreases with increasing distance from 
the center of the stellar disk. On the other hand, the intensity 
in the center of the line is significantly changed very close to 
the limb (see Fig. ll2i . These effects cannot be described by the 
simple formulae (|6jl and @. The difference between the line 
profile calculated by integrating the specific intensity over the 
disc and the profile which includes the velocity field and de- 
tailed radiative transfer is very small (see the magnified part of 
Fig.ll6i. as expected, since the radial velocity gradient is rela- 
tively small and the line forming region is very thin. For more 
extended and rapidly rotating sources these effects will be am- 
plified. 

3.4. Accretion discs 

Our method is also appropriate for solving the radiative transfer 
equation in accretion discs. It can take advantage of the accre- 
tion disc geometry. Since the disc is densest and, consequently, 
optically thickest close to the equatorial plane, we can employ 
the possibility of unevenly distributed angles 9 in the definition 
of the primary grid. 

In addition to having a better resolution of the disc, the cho- 
sen grid is also able to describe possible jets. Using this method 
we can also include in the calculations the boundary region, 
winds from this region and from the inner disc, as well as the 
hot corona beyond the disc. Another advantage of this method 
is the ability to handle both optically thin and thick discs. 
Emergent spectra from the accretion disk of a cataclysmic vari- 
able calculated using our method were presented in Korcakova 
et al. H2004b> . 



3.5. The tests of the grid 

The grid tests show a linear dependence of the computing time 
on the number of geometrical depth points D (see Fig. [n] 
left panel). The right panel of the same figure shows the de- 
pendence of the computing time on the number of frequency 
points N. The fitting function is a polynomial of the third or- 
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Fig. 17. The dependence of the computing time on the number 
of depth points {left panel) and on the number of frequency 
points (right panel). 
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Fig. 18. The dependence of the computing time on the number 
of angular grid points. 



der. Although one may expect only linear dependence of the 
computing time on a number of frequency points, higher order 
dependence is due to the necessity of interpolating the intensity 
for each frequency point at each boundary of the cells due to 
the Doppler shift. In Fig. we show the dependence of com- 
puting time on the number of angular grid points /. The fitting 
function is also a straight line. 

4. Conclusions 

We presented a new method for solving the radiative transfer 
equation in axial symmetry with the possibility of including ar- 
bitrary velocity fields. The basic idea is to solve the transfer 
equation in planes, that intersect the object. In a given plane a 
combination of long and short characteristics methods is used. 
This method allows us to better describe the global character of 
the radiation field and the resulting computing time is not too 
long. The velocity field is taken into account using the Lorentz 
transformation of frequency, which allows us to solve the trans- 
fer problem in the region with a small velocity gradient as well 
as for high velocity gradients. This technique is very useful for 
studying stellar wind, where it is applicable to the stellar wind 
region together with the stellar photosphere. 

Tests of this method were performed for a model atmo- 
sphere of a B type star with - 17 x 10 3 K, gravitational 
acceleration log g = 4.12, and radius 3.26R Q . For a static spher- 



D. Korcakova, J. Kubat: Radiative transfer in moving media II. 



11 



ically symmetric atmosphere our code gives correct results, as 
can be seen from Figure|9] We also present the limb darkening 
law for our model ('Fig.llOi. Note that the line profile shows 
limb brightening at the central line frequency (Fig. ll2> . This 
result is in agreement with the results of Hadrava and Kubat 
d2003l . 

Further tests were performed in the presence of a velocity 
field. For an expanding atmosphere (stellar wind) we adopt the 
classical /3- velocity law (|8}. The resulting line profile is shown 
in Fig.^] Since we take the input parameters (temperature, 
density) from the hydrostatic code, we do not obtain a P Cygni 
profile. However, the line profile is shifted to the blue part and 
deformed. Note that from these blue parts of ultraviolet lines 
one can determine the wind terminal velocity. 

In section l3"3l we studied the application of our method to 
the problem of stellar rotation. Fig. H6lshows the necessity of 
including the frequency dependence of limb darkening in cal- 
culating rotationally broadened line profiles. On the other hand, 
the difference between the flux calculated by integrating the 
Doppler shifted static profile across the disc and the "true" flux 
that included the solution of the transfer problem with a veloc- 
ity field is very small in a thin atmosphere. 

The dependence of computational time on the number of 
angular points (Fig.ll8> as well as on the number of geometrical 
points is linear (Fig. left panel). The dependence on the 
number of frequency points is more complicated, which stems 
from the interpolation of intensity due to the Doppler shift at 
cell boundaries. 

The limitation of this method is only in an axially symmet- 
rical approach and not very steep velocity gradient. In the latter 
case it is necessary to refine the grid and the computing time 
becomes longer. It is not possible to use it for the relativistic 
velocities too, since the aberration effect must be included in 
this case. 

Our method for solving the radiative transfer equation is es- 
pecially suitable for including stellar winds, since it is able to 
handle the outer region of a stellar wind together with the stellar 
photosphere. This is necessary, above all, for line-driven winds 
of hot stars. The process of initiating this type of wind near 
the photosphere is not fully understood yet. Calculations that 
involve the wind together with the quiet and possibly almost 
static stellar atmosphere may be able to resolve this problem in 
the future. Our method can also be used to solve the transfer 
problem in very rapidly rotating stars. Gravitational darkening 
is too large in these objects, so it is not possible to neglect the 
dependence of physical properties on the angle and to calculate 
using the assumption of spherical symmetry. The results of the 
application of our method can become useful for interpreting 
interferometric observations. The geometrical flexibility of our 
method is also very useful for studying extremely nonspherical 
objects such a accretion discs. It is possible to simultaneously 
include the disc, central object, jets and hot corona. Another 
advantage of this method is the possibility of solving both op- 
tically thin and thick discs. 
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